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INTRODUCTION: 

In  this  work,  we  are  seeking  to  construct  and  test  the  first  practical  full-field  transmission  ultrasound 
breast  imaging  system.  The  system  will  have  a  large  field-of-view  capable  of  producing  2D  images  of 
the  entire  breast  in  near  real-time.  Because  it  uses  a  fixed  imaging  geometry  that  does  not  involve 
scanning,  it  will  also  circumvent  the  operator-dependence  of  conventional  ultrasound  imaging  methods. 
The  compact  vertical  design  mimics  the  form  factor  of  mammography  systems  and  would  provide  a 
basis  for  later  incorporation  of  an  x-ray  source  and  detector,  which  would  allow  for  routine  dual¬ 
modality  x-ray  and  acoustic  imaging  for  screening  and  diagnosis.  If  successful,  this  work  will  provide 
the  basis  for  a  dual-modality  ultrasound  and  mammography  imaging  system  that  could  help  improve 
early  detection  of  breast  cancer,  especially  in  young  women  with  dense  breasts,  which  are  often  not  well 
visualized  mammographically.  It  would  also  help  reduce  false  positives,  which  add  to  the  expense  and 
anxiety  associated  with  current  approaches  to  breast  cancer  screening. 


BODY: 

Summary:  We  have  made  substantial  progress  in  this  first  year  of  the  project.  Most  tasks  in  the  original 
Statement  of  Work  for  year  1  have  been  completed  on  schedule.  The  only  delay  has  been  in  the 
optimization,  fabrication,  and  acceptance  testing  of  the  lenses,  which  we  had  scheduled  to  complete 
within  9  months  and  is  just  now  wrapping  up  in  the  13th  month.  This  arose  mainly  from  the  need  to  alter 
our  original  design  from  a  refraction-based  system  to  a  reflection-based  system,  due  to  unacceptably 
high  levels  of  acoustic  attenuation  predicted  in  the  refraction-based  system. 

Aim  1:  To  design  and  construct  a  high-resolution,  large  field  of  view  ultrasound  breast  imager  by 
combining  a  super  high-resolution  AO  sensor  and  a  large  aperture  acoustic  lens. 

Task  1:  Design,  fabricate,  and  test  sound  source 

l.a.  Simulate  acoustic  propagation  from  array  of  9  2.7”  elements  to  confirm  suitability  of  design  choice 
(U  of  C;  mos.  1-2). 

We  have  developed  theory  and  software  necessary  to  simulate  the  acoustic  field  emanating  from  square 
transducer  elements.  In  brief,  we  made  use  of  standard  Rayleigh-Summerfield  diffraction  theory  and  the 
angular  spectrum  decomposition  to  develop  expressions  for  the  field  propagating  from  a  rectangular 
piston-like  transducer.  By  the  linearity  of  the  resulting  equations,  we  can  use  superposition  to  determine 
the  field  produced  by  an  array  of  square  or  rectangular  transducers. 

The  details  are  provided  in  Appendix  A,  but  the  key  results  are  shown  in  Fig  1  below,  which  illustrates 
the  nonnalized  in-plane  acoustic  intensity  distribution  at  a  distance  2  inches  from  a  3”  by  3”  square 
transducer  piston  operating  at  a  CW  frequency  of  3.35  MHz.  The  figure  shows  both  a  surface  plot  and  a 
color  mapped  intensity  distribution. 


4 


I  4  5  0  5  4  « 


Figure  1:  (Left)  The  normalized  in-plane  acoustic  intensity  distribution  for  a  3”  by  3”  square 
transducer  piston  operating  at  a  CW  frequency  of  3.35  MHz.  The  pressure  field  intensity  is  depicted  in 
the  plane  parallel  to  the  transducer  surface  at  a  distance  of  2  ”,  and  a  background  medium  of  water  has 
been  assumed.  We  note  that  the  square  modulus  of  the  acoustic  pressure  field  is  plotted.  (Right)  This 
figure  presents  the  same  information  as  the  left  plot  but  shows  the  spatial  distribution  of  acoustic 
intensity  looking  at  the  plane  of  interest  head-on.  For  a  planar  detector  positioned  2  ”  from  the 
transducer  surface,  one  would  expect  to  observe  this  type  of  diffraction  pattern  for  an  image  acquired 
with  no  scattering  object  present. 


l.b.  Fabricate  source  array  (Santee;  mos.  2-3) 

A  large-area  sound  source  has  been  designed  and  fabricated  with  four  4  MHz,  3”x3”  size  elements.  The 
original  plan  called  for  working  at  3.35  MHz,  but  technical  difficulties  necessitated  a  change  to  4  MHz, 
which  should  not  affect  any  of  the  properties  of  the  planned  system.  A  quarter-wavelength  acoustic 
matching  layer  was  provided  to  allow  efficient  operation  in  water.  Each  element  was  electrically  tuned 
to  match  to  the  50-Ohm  impedance  of  an  RF  Amplifier  powered  by  a  4.0  MHz  electrical  signal  from  a 
digital  function  generator.  An  RF  power  splitter  was  also  designed  and  built  to  distribute  RF  power 
equally  to  all  4  elements  of  the  sound  source.  First-level  software  was  also  developed  to  allow  digital 
control  of  acoustic  power  output  of  the  sound  source.  The  current  design  is  being  extended  to  build  a 
larger  source  to  cover  the  8”x8”  field  area. 

1. c  Test  source  and  compare  with  simulation  results  (U  of  C;  mos.  3-4) 

Our  preliminary  comparisons  indicate  good  agreement  between  the  real  and  simulated  fields.  More 
detailed  comparison  work  is  ongoing. 

Task  2:  Fabricate  and  test  AO  sensor 

2.  a  Fabricate  1”  X  1”  AO  sensor  (Santee;  mos.  4-5) 

At  least  three  4  MHz  AO  sensors  have  been  built.  We  were  able  to  overcome  some  technical  challenges 
and  build  even  larger  sensors  than  originally  planned.  The  detection  area  of  each  sensor  was  at  least 
2”x2”.  The  acoustic  detection  sensitivity  of  these  sensors  was  found  to  be  in  the  milliWatt/cm2  range. 
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The  angular  acceptance  of  the  AO  sensor  is  at  least  20°  about  normal  incidence  condition.  The  acoustic 
detection  sensitivity  was  quite  dependent  on  the  acoustic  beam  angle,  and  might  impact  image  quality. 
However,  this  needs  to  be  studied  using  the  acoustic  lens  (or  mirrors)  being  designed  so  that  if  needed, 
the  AO  sensor  design  can  be  change  to  minimize  the  angular  effect.  The  milliWatt/cm2  range  acoustic 
detection  sensitivity  of  the  AO  sensor  should  be  adequate  to  allow  us  to  operate  the  Breast  Imager  below 
the  FDA  specified  safety  standard. 

2.b  Test  resolution  and  sensitivity  of  AO  sensor  (U  of  C;  mos.  5-6) 

We  sought  to  characterize  the  spatial  resolution  of  the  prototype  transmission  ultrasound  imaging  system 
by  deriving  estimates  for  the  edge  spread  function  (ESF),  the  line  spread  function  (LSF),  and  the  pre¬ 
sampled  modulation  transfer  function  (MTF)  of  the  system.  We  also  sought  to  characterize  the  noise 
properties  of  the  prototype  transmission  ultrasound  imaging  system  by  determining  an  estimate  for  the 
noise  power  spectrum  (NPS)  of  the  system. 

Spatial  Resolution  Analysis 

The  spatial  resolution  analysis  was  performed  using  the  angled  edge  technique  originally  proposed  by 
Reichenbach  et  al.  [1],  in  which  use  of  a  small  (~  3°)  edge  angle  pennits  finer  sampling  of  the  system 
ESF  than  the  pixel  pitch  of  a  digital  detector,  thus  minimizing  the  possibility  of  MTF  aliasing.  A 
transmission  ultrasound  image  of  a  slightly  angled  knife  edge  was  acquired  with  temporal  averaging  to 
reduce  noise.  A  segment  of  the  knife  edge  image  is  shown  in  Figure  2.  Each  row  of  pixels  in  Figure  2 
provided  a  shifted  sample  of  the  knife  edge.  The  individual  rows  were  spatially  registered  by  estimating 
the  edge  position  within 


Figure  2:  Image  of  the  knife  edge  used  to  compute  the  pre-sampled  MTF  of  the  imaging  system.  The 
image  was  acquired  using  a  peak  transducer  voltage  of  0. 7  V.  Continuous  frequency  sweeping  of  3.25- 
3.45  MHz  at  100  MHz/s  was  also  employed  to  minimize  speckle  artifacts  due  to  acoustic  wave 
coherence. 
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After  shifting  each  row  of  pixel  values  by  the  estimated  edge  location,  the  spatially  registered  data 
points  were  averaged  over  .0625-mm  intervals  to  obtain  an  estimate  of  the  system  ESF.  Following  the 
work  of  Boone  and  Seibert  [2],  the  measured  ESF  data  was  fit  to  a  combined  error  function  and 
exponential  fit  of  the  form 

ESF(x)  =  a  sgn(x  -  x0 ){ 1  -  exp(-  b|x  -  x0 1)}+  c  jerf  (Vd  (x  -  x0 ))}+  e. 


Analytic  estimates  for  the  normalized  system  LSF  and  MTF  were  determined  from  the  fit  parameters  a, 
b,  c,  d,  e,  and  xq  according  to  the  relations 


LSF(x)  = 


2(a  +  c) 


ab  exp(-  b|x  -  x0 1)+  2c.  —  exp(-  d(x  -  x0  )2 ) 


MTF  if) 


a  +  c 


cexp 


C  _2  rl\ 

n  J 

d 


+  a 


1  + 


4  ^2/2Vl 


Spatial  Resolution  Analysis  Results 

Figure  3  presents  the  measured  ESF  data  points  with  overlaid  exponential-error  function  regression 
model.  The  statistics  of  the  computed  fit  are  also  provided. 


Relative  Position  (millimeters) 


Figure  3:  Measured  ESF  data  and  analytic  regression  model  with  residual  subplot.  The  r2  value  for  the 
regression  curve  is  0.9999,  and  the  root-mean-square  error  (RMSE)  is  only  1.15%  of  the  mean 
measured  ESF  value.  These  statistics  indicate  that  the  combined  exponential-error  function  regression 
model  provides  an  exceiierit  fit  to  the  measured  data. 
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Figure  4  shows  the  normalized  analytic  LSF  and  the  analytic  pre-sampled  MTF  derived  from  the 
regression  equation  used  to  fit  the  ESF  data  in  Figure  6. 


Normalized  Analytic  LSF  for  Acoustomammography  System 


Analytic  Presampled  MTF  for  Acoustomammography  System 


Figure  4:  The  normalized  analytic  LSF  and  the  analytic  pre-sampled  MTF  for  the  transmission 
ultrasound  imaging  system.  These  curves  were  obtained  from  the  fitting  parameters  of  the  combined 
exponential- err  or  function  regression  model  used  to  fit  the  measured  ESF  data.  The  full-width-at-half- 
maximum  (FWHM)  of  the  system  LSF  indicates  an  approximate  spatial  resolution  of 400  microns. 

The  FWFIM  of  the  analytic  LSF  was  calculated  to  be  0.3997  millimeters,  indicating  an  approximate 
spatial  resolution  of  400  microns  for  the  prototype  transmission  ultrasound  imaging  system.  The  analytic 
curve  for  the  pre-sampled  MTF  shows  a  typical  Gaussian  form.  The  curve  is  equal  to  unity  at  zero 
spatial  frequency  and  falls  off  smoothly  to  nearly  zero  at  the  system  Nyquist  frequency  (3.79 
cycles/mm). 

Noise  Analysis 

The  NPS  of  the  imaging  system  was  determined  by  acquiring  25  detector  flood  images  using  the  same 
transducer  settings  for  each  acquisition.  One  such  image  is  included  in  Figure  5.  For  each  flood  image,  a 
subsection  of  the  total  image  matrix  was  extracted  over  which  the  detector  illumination  appeared  highly 
unifonn.  The  NPS  of  each  subsection  was  computed  by  subtracting  the  mean  pixel  value  <PV>  and 
computing  the  square  modulus  of  the  two-dimensional  discrete  Fourier  transfonn  of  the  resulting  matrix. 
In  the  continuum, 


NPS(f) 


j  j[PV(r)-<PV  >\e~j27ctrd: 


The  25  separately  computed  noise  power  spectra  were  averaged  to  obtain  a  mean  estimate  for  the  NPS 
of  the  system.  This  estimate  was  nonnalized  by  its  largest  amplitude  component  and  plotted  as  a  two- 
dimensional  function  of  spatial  frequency  f. 
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Figure  5:  Example  of  one  of  the  detector  flood  images  used  to  estimate  the  NPS  of  the  imaging  system. 
The  image  shows  highly  uniform  illumination  of  the  AO  detector  field  of  view.  A  peak  transducer  voltage 
of  2.6  V  and  frequency  sweeping  of 3.25-3.45  MHz  at  a  rate  of  100  MHz/s  were  used  because  of  the 
uniform  detector  illumination  provided  by  these  settings 

Figure  6  presents  an  estimate  of  the  two-dimensional  NPS  of  the  imaging  system  as  a  surface  plot.  The 
height  of  the  NPS  indicates  the  noise  power  contained  per  unit  frequency  bandwidth  at  each 
corresponding  spatial  frequency.  It  can  be  seen  that  significant  noise  power  is  contained  at  lower  spatial 
frequencies  near  zero  frequency.  The  NPS  amplitude  falls  steadily  when  moving  toward  higher  spatial 
frequencies. 
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Verification  of  "1/f"  Noise  Power  Spectrum  for  Acoustomammography  System 


Figure  6:  Surface  plot  of  the  two-dimensional  NPS  of  the  transmission  ultrasound  imaging  system.  The 
lower  two  graphs  show  single  slices  taken  through  the  NPS  surface  plot  along  the  x-  and  y-spatial 
frequency  axes.  These  slices  are  plotted  on  log-log  axes  to  demonstrate  the  presence  of  1/f  noise  in  the 
imaging  system. 


The  lower  two  plots  in  Figure  6  show  single  slices  of  the  NPS  taken  along  the  x-spatial  frequency  and  y- 
spatial  frequency  axes.  For  these  graphs,  the  NPS  is  plotted  on  log-log  axes.  These  slices  are  well  fit  by 
power  law  regression  models,  which  appear  as  straight  lines  on  the  log-log  plots.  The  characteristic 
exponent  values  of  the  regression  models  are  0.680  mm’1  and  0.665  mm’1  for  the  x-  and  y-directions, 
respectively.  This  type  of  behavior  is  ubiquitous  in  nature  and  is  termed  “1  If  noise.”  It  is  caused  by 
random  fluctuations  in  the  solid-state  circuit  elements  comprising  most  electronic  devices. 

Task  conclusions 

•  Though  the  intrinsic  spatial  resolution  of  the  AO  detector  is  itself  very  high  (~  10’9  meters)  the 
spatial  resolution  of  the  imaging  system  is  substantially  reduced  by  the  pixelation  of  the  digital 
camera  used  to  record  the  display  of  the  AO  detector  for  image  processing  and  storage.  The 
measured  resolution  is  400  microns. 
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The  noise  present  in  the  imaging  system  is  primarily  electrical  in  nature,  as  evidenced  by  the  1  If 
noise  behavior  observed  in  the  system  NPS  in  both  the  x-  and  y-spatial  directions. 


Task  3:  Incorporate  video  camera 

3.  a.  Incorporate  the  new  video  camera  with  source  and  AO  sensor  (no  lens  yet)  and  test  resolution  and 
sensitivity  (U  of  C;  mo.  6) 

We  have  not  yet  incorporated  the  new  video  camera.  The  results  above  were  acquired  with  the  original 
video  camera.  Adding  the  video  camera  will  be  a  simple  modification  to  incorporate  along  with  the 
lenses  but  the  optimal  specifications  needed  for  the  camera  will  depend  on  the  properties  of  the 
fabricated  lenses. 

Task  4:  Design,  fabricate,  and  test  acoustic  lens 

4.  a.  Determine  final  lens  specifications  based  on  simulation  studies  of  whole  system  (U  of  C;  mos  1-3) 
Our  initial  plan  was  to  fabricate  an  afocal  triplet  of  acoustic  lenses,  as  shown  in  Fig.  7. 


Breast 

side 


Low  F /#  is 
required 


Field  lens 


Equivalent 

Petzval 

Type 

design 


Fig.  7:  Initial  afocal  triplet  design.  The  acoustic  attenuation  in  this  design  would  have  been  too  high. 


Detector 

side 


However,  initial  design  and  simulation  studies  showed  that  because  of  the  fast  speed  of  the  lens  required 
(~F/3  or  less),  the  lens  thickness  variation  would  have  been  about  ~60  mm  which  would  have  caused 
unacceptable  attenuation  of  the  acoustic  field. 


11 


An  alternative  single-lens  design,  shown  in  Fig.  8,  would  still  have  been  too  thick  (56  mm),  leading  to 
two  or  three  orders  of  magnitude  of  acoustic  attenuation.  Therefore  it  was  decided  to  produce  a 
reflective  system  rather  than  a  refractive  lens. 


Fig  8:  A  single-lens  design  aimed  at  reducing  acoustic  attenuation  still  would  have  had  unacceptable 
levels  of  it.  We  turned  instead  to  a  reflection-based  system. 

4.b.  Optimize  lens  design  using  Zemax  lens  design  software  (U  of  C/Sasian;  mos.  4-6) 

We  used  the  Zemax  lens  design  software  to  optimize  the  mirror-based  system.  The  schematic  of  the 
optimized  mirror-based  system  is  shown  below  in  Fig  9. 

Secondary 

mirror 


Small  attenuation  ~0.2 

Fig  9:  Cross-sectional  schematic  of  the  lens-based  system. 


Detector 

side 
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A  three-dimensional  rendering  is  shown  in  Fig.  10. 


Fig  10:  Three-dimensional  view  of  the  mirror-based  system. 
The  system  characteristics  are  as  follows: 

•  Length:  600  mm 

•  Width:  310  mm 

•  Working  F/#  =  2.88 

•  Wavelength:  0.5  mm 

•  Stop  aperture  at  primary  mirror 

•  Other  apertures  for  stray  sound  suppression  are  allowed 

•  System  is  plane  symmetric  in  y-z  plane 

•  Magnification  0.5 

The  primary  mirror  surface,  which  is  concave,  is  described  by 


sag 


with  parameters 

•  c  = -1/1000  mm 

•  A  =  1 .55586 1  e-4  mm'1 

•  B  =  5.74993e-7  mm'2 

•  C  =  5.160419e-7  mm'2 

•  Diameter  =130  mm 


With  these  specifications,  the  expected  performance  is  given  in  the  Airy  disk  diagrams  shown  in  Fig.  1 1 . 
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Fig  11:  Airy  spot  diagrams  showing  expected  performance  of  the  lens  system. 


4.c.  Fabricate  lens  (U  of  C/Sasian;  mos.  7-8) 


The  mirrors  are  being  fabricated  by  Lam  Optics  in  Tucson,  AZ.  The  secondary  mirror  is  finished,  as  is 
the  first  toroidal  primary.  The  second  primary  mirror  is  expected  to  be  finished  by  June  30th. 

4.d.  Perform  acceptance  testing  on  lens  (U  of  C;  mo.  9) 

Acceptance  testing  will  follow  final  delivery  of  the  lens,  but  our  consultant  Jose  Sasian  has  been 
personally  visiting  the  fabrication  facility  to  observe  fabrication  and  oversee  in-house  quality  assurance. 


Task  5:  Assemble  complete  system 

5. a.  Fabricate  immersion  tank  (Santee;  mo.  9) 

The  dimensions  of  the  immersion  tank  will  be  determined  by  the  final  configuration  of  lenses,  which 
must  await  final  acceptance  testing  of  the  lenses. 

5.b.  Assemble  system  (U  of  C;  mo.  10-12) 

Assembly  must  obviously  await  the  delivery  of  the  lenses. 
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KEY  RESEARCH  ACCOMPLISHMENTS: 


•  Successful  simulation  and  construction  of  acoustic  field  of  multielement  sound  source. 

•  Successful  construction  of  2-inch  acousto  optic  (AO)  sensor. 

•  Measurement  of  noise  and  spatial  resolution  properties  of  AO  sensor  coupled  with  existing  video 
camera. 

•  Development  and  characterization  of  multiple  acoustic  lens  systems,  resulting  in  a  very  original 
acoustic-mirror  design  that  should  minimize  acoustic  attenuation  relative  to  the  original 
refraction-based  designs. 

REPORTABLE  OUTCOMES: 

•  Poster  presentation:  J.R.  Rosenfield,  J.S.  Sandhu,  J.K.  Tawiah,  andP.J.  La  Riviere,  “Evaluation 
of  the  Spatial  Resolution  and  Noise  Properties  of  a  Prototype  Transmission  Ultrasound  Breast 
Imaging  System  Employing  an  Acousto-Optic  Detector,”  NIH  NIB  IB  training  grant  symposium, 
June,  2012. 

•  Oral  presentation:  J.R.  Rosenfield,  J.S.  Sandhu,  J.K.  Tawiah,  andP.J.  La  Riviere,  “Evaluation  of 
the  Spatial  Resolution  and  Noise  Properties  of  a  Prototype  Transmission  Ultrasound  Imaging 
System  Employing  An  Acousto-Optic  (AO)  Detector,”  AAPM  Annual  Meeting,  July,  2012. 


CONCLUSION: 

In  summary,  we  have  made  excellent  progress  through  the  year.  The  simulation  and  construction  of  the 
sound  sources  and  AO  sensors  have  gone  according  to  schedule  and  careful  measurements  have  shown 
that  the  system  resolution  of  the  AO  sensor  coupled  with  the  current  video  sensor  is  -400  microns.  The 
noise  present  in  the  imaging  system  is  primarily  electrical  in  nature,  as  evidenced  by  the  1  If  noise 
behavior  observed  in  the  system  NPS  in  both  the  x-  and  v-spatial  directions.  Our  initial  lens  design 
proved  impractical  because  of  the  large  acoustic  attenuation  of  the  materials  needed  to  achieve  the 
design,  but  we  have  produced  a  strong  alternative  design  based  on  the  use  of  acoustic  mirrors. 
Fabrication  is  going  well  and  should  be  complete  this  summer.  Integration  and  testing  of  the  complete 
system  will  then  proceed  as  originally  planned. 

REFERENCES: 

[1]  S.  E.  Reichenbach,  S.  K.  Park,  and  R.  Narayanswamy,  "Characterizing  digital  image  acquisition 
devices,"  Opt.  Eng.  (Bellingham)  30,  170-177  (1991). 

[2]  J.M.  Boone  and  J.  A.  Seibert,  “An  analytical  edge  spread  function  model  for  computer  fitting  and 
subsequent  calculation  of  the  LSF  and  MTF,”  Med  Phys.  21,  1541-1545  (1994). 
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APPENDICES: 


Appendix  A  provides  significant  detail  on  the  acoustic  simulation. 

Appendix  A:  Details  of  acoustic  simulation 


I.  Wave  Equation  Description 

The  imaging  system  relevant  to  ultrasound  diffraction  tomography  can  be  treated  as  an 
inhomogeneous  medium  of  finite  size  ( i.e .,  the  imaged  object)  immersed  within  a  lossless,  homogeneous 
background  medium  providing  acoustic  coupling.  A  schematic  of  the  measurement  geometry  is  given  in 
Figure  1.  In  the  calculations  that  follow,  all  waves  considered  will  be  compressional  {i.e.,  longitudinal) 
in  nature,  as  opposed  to  shear  waves  that  may  be  present  in  viscous  materials.  Furthermore,  for  the  sake 
of  generality,  attenuation  of  mechanical  ultrasound  energy  within  the  imaged  object  will  be  implicitly 
assumed.  The  fonnalism  developed  here  will  be  general  enough  that  the  case  of  pure  scattering  of  an 
incident  wave  field  by  the  object  with  no  a  ttenuation  could  just  as  well  be  described  by  the  same 
equations. 

For  inhomogeneous  media,  the  propagation  of  the  acoustic  pressure  field  in  a  source-free  region 
is  well  known  to  satisfy  the  non-dispersive  wave  equation 


vw, 

c  (r )  ot 


(1.1) 


where  u(f  ,t)  represents  the  complex  amplitude  of  the  field  as  a  function  of  position  r  and  time  t,  and 
c(r)  denotes  the  velocity  of  sound  as  a  function  of  position.  In  an  inviscid  medium,  the  wave  velocity 


c(r)  is  a  purely  real  quantity,  but  in  the  presence  of  attenuation,  the  velocity  becomes  complex  with  the 


imaginary  part  reflecting  attenuation  in  the  supporting  medium.  The  operator  V2  denotes  the  Laplacian 
for  the  spatial  coordinates  of  the  field  point  r  .  In  order  to  determine  the  time-  and  space-dependent 
acoustic  pressure  field  produced  by  a  vibrating  source,  we  must  solve  (1.1)  for  the  manner  in  which  the 
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acoustic  field  is  propagated,  and  we  must  constrain  the  solution  to  match  the  initial  or  boundary 
conditions  defined  by  the  problem. 
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Figure  1:  The  measurement  geometry  of  ultrasound  diffraction  tomography.  The  relative  positions  of  the  source  and  detector 
planes,  the  transducer  aperture,  and  the  imaged  object  are  shown.  The  object  is  immersed  in  the  lossless  background  medium. 
All  waves  considered  will  be  compressional  (i.e.,  longitudinal)  in  nature,  as  opposed  to  shear  waves  that  may  be  present  in 
viscous  materials.  For  the  sake  of  generality,  attenuation  of  mechanical  ultrasound  energy  within  the  imaged  object  will  be 
assumed. 

Taking  the  Fourier  transform  of  (1.1)  with  respect  to  time,  we  find  that  the  equation  conjugate  to  (1.1)  is 

V2u(f,t)  +  k2(f)u(f,t)  =  0,  (1.2) 

where  k(f )  is  the  wave  number  as  a  function  of  position  for  the  acoustic  field  having  angular  frequency 
CO .  The  wave  number  k(r  )  is  in  general  a  complex  quantity,  with  the  imaginary  part  representing  any 
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attenuation  of  ultrasound  energy  that  may  be  present  in  a  supporting  medium  and  the  real  part  specifying 
the  phase  velocity  of  the  field  at  a  particular  position. 


In  the  lossless  background  medium,  lc(f )  is  a  real  number  having  the  magnitude 


^  _  2/r  _  co 


(1-3) 


where  \  and  c0  are  the  wavelength  and  velocity  of  the  field  in  the  background  medium,  respectively. 

We  note  that  for  our  application,  we  will  only  consider  a  continuous  wave  (CW)  acoustic  source 
operating  at  a  single  vibrational  frequency,  so  that  CO  corresponds  to  the  frequency  of  the  transducer 
used. 

As  (1.2)  is  expressed  in  terms  of  a  single  angular  frequency  CO,  we  may  write  the  functional 
dependence  of  the  acoustic  pressure  phasor  u(r,t )  as  the  product  of  a  spatially  varying  function  and  a 
time-harmonic  factor  thus: 

u(r  ,t)  =  u{r)eJC0t .  (1.4) 


Upon  substitution  of  this  expression  into  ( 1 .2),  it  can  immediately  be  seen  that  the  time  dependence  of 
the  field  is  suppressed,  and  one  obtains  for  the  spatial  dependence 

\/2u(r)  +  k2(r)u(r)  =  0.  (1.5) 

At  this  point,  it  is  customary  to  write  the  wave  number  k(r )  in  terms  of  the  complex  acoustic  refractive 
index  n(r )  of  the  inhomogeneous  medium.  In  particular,  if  we  denote  by  ng(r )  the  refractive  index 
deviations  of  the  scattering  object  from  the  homogeneous  background  medium,  then  we  have 

k{r)  =  k0n(r )  =  k0  (1  +  ns  (r )),  (1.6) 

where  k0  represents  the  wave  vector  magnitude  in  the  background  medium.  It  should  be  emphasized 
that  (1.6)  is  valid  both  within  and  outside  of  the  scattering  object.  For  locations  outside  the  scattering 
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object,  ns(r )  =  0  and  (1.6)  reduces  to  k (r )  =  k0 .  For  completeness,  we  note  that  the  acoustic 
refractive  index  may  be  written  as 


n(f)  = 


c(r) 


(1.7) 


where  C0  again  denotes  the  speed  of  sound  in  the  homogenous  background  medium  in  which  the 
imaged  object  is  immersed. 

To  put  (1.5)  into  a  more  enlightening  form,  we  observe  that  if  we  substitute  (1.6)  into  (1.5),  add  a 
factor  of  k^u(r  )  to  both  sides,  and  rearrange,  we  obtain 


{V-  +kl)u{r)  =  -k;{n\r)-\)u{r). 


(1.8) 


If  we  consider  that  the  acoustic  field  at  r  consists  of  an  incident  field  w0(r)  present  without  the 
scattering  object  in  place  (i.e.,  the  field  due  solely  to  the  output  of  the  acoustic  transducer)  and  a 
scattered  field  contribution  Us(r )  due  entirely  to  the  presence  of  the  inhomogeneous  scattering  object, 
then  we  may  write 

u(r)  =  u0(r)  +  us(r ).  ( 1 .9) 

This  division  of  the  total  field  u(r )  into  an  incident  component  and  a  scattered  component  allows  us  to 
write  (1.8)  as  the  following  equivalent  system  of  equations: 

(V2  +  kl)uQ{r)  =  0. 

<  (V2  +kl)us{r)  =  ~u{r)o{r).  (1.10) 

u(r)  =  u0(r)  +  us(r). 

Here,  o(r )  denotes  the  complex  scattering  object  function,  defined  as 

o{r)  =  kl{n\r)-\).  (1.11) 

Note  that  o(f )  is  defined  both  within  and  outside  of  the  scattering  object,  although  it  vanishes 

everywhere  outside  of  the  scattering  object  since  according  to  (1.6),  n(r )  =  1  in  the  homogeneous 
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background  medium.  A  final  technical  point  worth  noting  is  that  the  acoustic  refractive  index  function 
n(f )  is  in  general  a  complex -valued  quantity.  For  lossless  media,  n(r)  is  real.  The  complex  nature  of 

the  refractive  index  for  attenuating  media  plays  an  extremely  important  role  in  ultrasound  diffraction 
tomography  using  measurements  of  field  intensity. 

II.  The  Angular  Spectrum  Method 

The  physical  significance  of  w0(r)  can  be  particularly  well  understood  by  noting  that  the 
differential  equation  in  (1.10)  for  the  incident  field  M0(r)  , 

(V2+£0>0(r)  =  0,  (1.12) 

is  the  homogeneous  Helmholtz  equation  for  the  acoustic  pressure  field  in  a  uniform,  lossless  medium 
characterized  by  a  wave  vector  of  magnitude  k0 .  The  solution  to  this  equation  is  simply  the  solution  to 

the  wave  equation  in  the  absence  of  the  scattering  object,  which  is  what  we  intended  W0(r )  to  represent 
when  we  defined  the  total  field  u(r)  according  to  (1.9).  The  solution  to  (1.12)  for  the  incident  field 
M0(r)  can  easily  be  verified  through  substitution  to  be  a  plane  wave 

u0(r)  =  efi*  (U3) 

moving  in  an  arbitrary  direction  specified  by  the  propagation  vector 

k  =(kx,ky,kz),  (1.14) 

which  has  the  magnitude  k  —  ^k]  +  k~  +  k 2  =  k0  appropriate  for  the  background  medium. 

Equation  (1.13),  together  with  the  condition  k  =  J k J  +  k2y  +  k2  =  k() ,  provides  the  generic 

form  of  the  plane  wave  solution  to  the  homogeneous  Helmholtz  equation  (1.12)  for  the  incident  acoustic 
field.  Because  (1.12)  is  a  linear  differential  equation,  a  superposition  of  plane  wave  solutions  each 
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having  the  form  given  in  (1.13)  and  subject  to  k  =  +  ky  +  k~  —  k0  is  therefore  also  a  solution.  In 

other  words,  we  can  construct  an  arbitrary  wave  field  that  satisfies  (1.12)  from  a  sum  of  plane  waves  all 
having  the  same  wave  vector  magnitude  k0  but  propagating  in  different  directions.  This  observation 
forms  the  mathematical  basis  of  the  so-called  angular  spectrum  method  (ASM)  for  solving  the 
homogeneous  Helmholtz  equation  subject  to  a  particular  set  of  boundary  conditions. 

In  three  dimensions,  the  ASM  enables  the  incident  acoustic  pressure  field  to  be  determined 
throughout  all  space  from  knowledge  of  the  field  on  only  a  single  plane.  To  understand  how  this  is 
possible,  we  consider  our  ultrasound  tomography  system  to  consist  of  a  square  piston  transducer  plate  of 
dimensions  \a,a\  centered  on  the  z-axis  and  lying  entirely  in  the  x-y  plane.  We  suppose  that  this 
transducer  is  producing  an  acoustic  wave  field  only  in  the  positive  z-direction.  In  the  plane  z  =  0  (/.  e.,  on 
the  transducer  surface),  the  incident  acoustic  pressure  field  U0(r )  can  be  decomposed  via  the  non¬ 
unitary,  two-dimensional  inverse  Fourier  transform  relation 

1  OO  CO 

u0(x,y,z  =  0)  =  — —  J  J a{kx,ky : z  =  0 )eJ(Kx+ky>)dkxdky,  (1.15) 

■  ' *  —00—00 

where  a( kx,k  :z  =  0)  denotes  the  spatial  frequency  spectral  density  function  required  to  synthesize 
the  known  pressure  distribution  u0(x,y,z  =  0)  from  two-dimensional  plane  waves  of  the  form 

j i^k  x ~\~k  y j 

e  '  y  .  Equation  (1.15)  expresses  the  so-called  angular  spectrum  decomposition  of  the  incident 
pressure  field  U0(r )  on  the  plane  z  =  0.  That  such  a  decomposition  can  be  achieved  is  a  fundamental 

result  of  Fourier  analysis.  Provided  U0(x,y,z  —  0)  is  known  for  all  x  and  y  in  the  plane  z  =  0,  we  can 
detennine  the  precise  functional  form  of  cc{kx,k  z  =  0)  required  to  synthesize  this  spatial  pressure 
distribution  simply  by  inverting  (1.15)  thus: 

co  oo 

a(kx,ky :  z  =  0)  =  J  J u0(x,y,z  =  0)e  J(kiX+ky} ]dxdy.  (1.16) 
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We  note  that  each  two-dimensional  plane  wave  present  in  the  summation  of  (1.15)  can  be  attributed  to 
one  of  the  valid  plane  wave  solutions  (1.13)  to  the  homogeneous  Helmholtz  wave  equation  (1.12).  In 
particular,  we  observed  earlier  that  any  plane  wave  of  the  form 


and  subject  to  the  condition 


u0(r)  =  e 


j(kxx+kyy+kxz ) 


(1.17) 


K  =  K  +  K  +  K 


(1.18) 


j  ^  x~\~k  j 

satisfies  the  differential  equation  (1.12).  Thus,  we  see  that  each  plane  wave  of  the  fonn  e  y 
appearing  in  the  summation  of  (1.15)  can  be  associated  with  a  valid  plane  wave  solution  eJ(k'x+  y>+  zZ>  to 
the  homogeneous  Helmholtz  equation,  where  the  value  of  Ay  is  determined  uniquely  from  the  values  of 
kx  and  k  according  to 


(1.19) 


In  (1.19),  we  have  restricted  values  of  Ay  to  be  positive  because  the  transducer  for  our  system  produces 
an  acoustic  field  only  in  the  positive  z-direction. 

Once  the  decomposition  of  the  incident  acoustic  field  in  the  transducer  plane  (z  =  0)  has  been 
achieved  via  (1.16),  we  can  detennine  the  field  in  any  other  parallel  plane  z  >  0  in  the  following  manner. 

Each  plane  wave  ej(k'x+k>’y+  zZ)  required  at  z  =  0  to  synthesize  U0(x,y,z  =  0)  according  to  (1.15)  will 

undergo  a  phase  shift  as  it  propagates  from  the  transducer  plane  to  the  plane  of  interest.  The  change  in 
phase  in  passing  from  the  point  (x,  y,  0)  in  the  transducer  plane  to  the  point  (x,  y,  z)  in  the  measurement 

plane  is  simply  Ayz.  Thus,  the  complex  amplitude  of  the  plane  wave  component  e'<k'x+  >y)  of  the  total 


acoustic  field  at  z  =  0  is  related  to  its  complex  amplitude  at  z  >  0  by  the  factor  ejk~~z .  This  is  just  another 
way  of  saying  that  for  a  given  plane  wave  component  of  the  total  acoustic  field, 

«,(i,;,z)  =  ac'(wl  =(ae',w)c"  =u0(x,y,  0)eA’.  (1.20) 

23 


Accounting  for  the  change  in  phase  of  each  plane  wave  component  of  the  angular  spectrum 
decomposition  of  U0(x,y,z  =  0)  in  propagating  from  the  plane  z  =  0  to  any  plane  z  >  0,  we  therefore 

have  for  the  acoustic  pressure  field  u0(x,y,z) 


w0 (x, y, z )  =  —  J  J a(kx , ky  :z  =  0)e 


](kxx+kvy+kzz) 


dkxdky , 


(1.21) 


where  again 


(1.22) 


as  required  by  the  homogeneous  Helmholtz  equation  for  the  background  medium.  Equation  (1.21) 
allows  us  to  calculate  the  incident  acoustic  field  at  any  point  (x,  y,  z)  by  detennining  a(kx,k  z  =  0)  , 


which  itself  requires  knowledge  only  of  U0(x,y,z  =  0)  .  We  note  that  we  can  write  (1.21)  in  terms  of 
the  spatial  frequency  spectral  density  function  cc{kx,  k  z)  for  any  plane  of  constant  z,  not  necessarily 
z  =  0,  by  observing  that 

a(kx ,  ky :  z)  =  a(kx ,  ky :  z  =  0)eJkzZ .  ( 1 .23) 

Then  (1.21)  assumes  the  recognizable  form  of  a  simple  two-dimensional  inverse  Fourier  transfonn  in  the 
plane  z: 

1  OO  00 

«0  (w  y> z )  =  z  f  \  j  a(kx » K  :  z)eJ(kxX+k,y)dkxdky .  ( 1 .24) 

—CO—OO 

Before  applying  the  ASM  to  a  rectangular  piston  transducer  having  a  known  surface  velocity 
profile,  it  is  important  to  consider  two  limiting  cases  for  the  z-component  kx  of  the  propagation  vector 

k  —  ( kx,k  ,k. )  .  We  mentioned  above  that  for  each  plane  wave  in  the  expansion  (1.21),  the  value  of  kz 
for  that  component  of  the  expansion  is  determined  from  the  corresponding  values  of  kx  and  k , 
according  to  (1.22).  Because  the  integration  limits  in  (1.21)  run  from  —  oo  to  +  oo ,  we  see  that  for 
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values  of  k  v  and  ky  such  that  <  k  +  kv ,  the  radical  in  (1.22)  becomes  complex  and  kz  is  pure 
imaginary.  Under  this  circumstance,  the  plane  wave  solution  to  the  Helmholtz  equation  becomes  the  so- 
called  evanescent  wave.  The  complex  nature  of  k _  for  evanescent  waves  leads  to  areal,  attenuating 

component  that  causes  the  amplitudes  of  these  solutions  to  decay  exponentially.  Though  they  are  valid 
solutions  to  the  wave  equation,  evanescent  waves  carry  no  ne  t  energy  and  can  be  considered  a  non¬ 
physical  outcome  of  the  ASM  related  to  boundary  condition  matching.  Evanescent  waves  decay  rapidly 
when  moving  away  from  any  boundaries  present,  and  their  influence  can  often  be  ignored  at  distances 

from  boundaries  exceeding  ten  times  the  acoustic  wavelength.  On  the  other  hand,  for  k°~  >  k]  +  k1 ,  k 

is  real  and  the  corresponding  solution  to  the  Helmholtz  equation  remains  a  plane  wave  propagating  in 
the  direction  of  positive  z  without  exhibiting  evanescent  decay. 

III.  Solution  for  the  Incident  Field  uo(r) 

We  will  now  apply  the  ASM  outlined  above  to  our  square  transducer  plate  having  a  known 
velocity  profde.  For  a  planar  acoustic  source  radiating  into  an  unbounded,  lossless  background  medium, 
solving  the  Helmholtz  wave  equation  (1.12)  for  the  incident  field  requires  matching  the  boundary 
condition  for  the  pressure  field  on  the  acoustic  source  plane  only.  Since  there  are  no  other  physical 
boundaries  present,  the  value  of  the  field  at  all  other  points  in  the  homogeneous  medium  is  detennined 
by  the  movement  of  the  wave  field  output  from  the  source  plane  through  the  medium  as  required  by  the 
propagation  equation  (1.1).  One  could  nominally  require  that  the  magnitude  of  the  acoustic  pressure 
field  vanishes  at  an  infinite  distance  from  the  transducer  surface.  We  will  see,  however,  that  this 
condition  will  be  satisfied  by  the  solution  we  obtain  by  only  considering  the  boundary  condition  on  the 
transducer  plane.  We  discussed  in  the  previous  section  how  the  ASM  enables  us  to  achieve  this  required 
boundary  condition  on  the  source  plane,  and  that  the  total  acoustic  field  can  be  determined  on  any  other 
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parallel  plane  simply  by  propagating  the  plane  wave  solutions  of  the  Helmholtz  equation  from  the 
source  plane  to  the  plane  of  interest. 

As  mentioned  previously,  our  analysis  will  be  confined  to  an  acoustic  source  consisting  of  a 
single  square  transducer  plate  of  dimensions  \a,  a\  centered  on  the  z-axis  and  lying  entirely  in  the  x-y 
plane.  We  will  assume  CW  excitation  of  the  source  at  a  single  vibrational  frequency  CO .  The  source  will 
be  embedded  in  an  ideal  rigid  baffle  such  that  there  is  no  vibration  of  the  transducer  plane  outside  of  the 
transducer  aperture  dimensions  \a,  Cl\.  We  can  therefore  write  for  the  surface  velocity  of  any  point  in 
the  transducer  plane 


v:(x,y,z 


0  :t)  =  v0rect 


rect 


(V| 
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-jcot 


(1.25) 


where  V0  is  a  constant  indicating  the  velocity  amplitude  of  the  oscillations  and  the  subscript  z  reminds 

us  that  this  is  the  velocity  component  perpendicular  to  the  plane  z  =  0.  Equation  (1.25)  is  known  as 
the  vibration  velocity  wavefonn  of  the  acoustic  source  aperture,  and  in  particular,  the  function 


g(x,y)  =  rect 


rect 


\aj 


fZl 

\a) 


(1.26) 


describing  the  spatial  dependence  of  the  velocity  waveform  is  known  as  the  apodization  function.  For 
clarity,  it  should  be  realized  that  for  a  moveable  surface  enclosing  a  fluid  acoustic  medium,  kinematic 
boundary  conditions  require  that  the  component  of  the  fluid  velocity  normal  to  the  surface  equals  the 
velocity  of  the  bounding  surface  in  the  same  direction.  In  other  words,  (1.25)  equally  well  represents 
either  the  velocity  of  the  vibrating  transducer  surface  or  the  component  of  the  fluid  velocity  normal  to 
the  transducer  surface  in  the  plane  z  =  0. 

Solution  of  the  Helmholtz  equation  (1.12)  requires  that  either  the  acoustic  pressure  on  t  he 
boundary  surface  So  or  the  normal  component  of  the  medium  particle  velocity  on  the  boundary  surface 
be  specified.  Specification  of  the  nonnal  component  of  the  medium  particle  velocity  distribution  in  the 
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transducer  plane  constitutes  the  so-called  Neumann  boundary  condition  for  the  differential  equation 
(1.12),  whereas  specification  of  the  pressure  distribution  on  this  plane  constitutes  the  Dirichlet  boundary 
condition.  As  written,  ( 1 .24)  is  expressed  in  terms  of  the  pressure  distribution  in  the  transducer  plane 
and  is  therefore  related  to  the  Dirichlet  boundary  condition.  However,  the  Neumann  boundary  condition 
is  more  directly  related  to  the  velocity  wavefonn  of  the  transducer  aperture,  and  therefore  it  has  proven 
advantageous  to  express  ( 1 .24)  in  terms  of  the  Neumann  condition.  This  can  be  achieved  by  first  taking 
the  two-dimensional  Fourier  transform  of  (1.24)  with  respect  to  x  and  y  to  obtain 

co  oo 

a(kx,ky  :z)=  J  J u0(x,y,z)e  ,lkrX+kf'v>dxdy.  (1.27) 

-00-00 

Equation  (1.27)  can  equally  well  be  expressed  in  terms  of  the  velocity  potential  function  ®0(r,?)  for 
the  incident  field,  which  for  an  inviscid  medium  is  related  to  the  incident  acoustic  pressure  field 
u0(r,t )  according  to 

d®0(7, 

dt 

Here,  p0  is  the  mass  density  of  the  uniform,  lossless  background  medium.  With  u0(r,  t )  =  u0(f)e 
as  in  ( 1 .4),  we  therefore  have 

®0(r,0  =  —  u,{r)eJ~  =®0(r)e~J*-  (1-29) 

JWPo 

Canceling  the  harmonic  factor  e  JC0>  from  both  sides  of  (1.29)  that  characterizes  the  time  dependence  of 
®0(r,?)  and  u0(f,t),  we  have 


u0(r,t)  =  p0 


t ) 


(1.28) 


- - «0(r). 

ja>Po 


Substitution  of  ( 1 .30)  for  u0  (r )  in  ( 1 .27)  gives 


(1.30) 


27 


a(kx,ky :  z)  =  -ja>p0  J  J 0(x,y,z)e  AKx+kyy) dxdy.  (1.31) 


Differentiation  of  ( 1 .3 1)  with  respect  to  z  gives 


da(kx,k  :z)  .  }}dOQ(x,y,z)  } 

- — - =  -JG>Po\\ - - - e  ”  dxdy,  (1.32) 


while  differentiation  of  (1.23)  with  respect  to  z  similarly  gives 

da(k  ,k  :z)  d(ejkzZ )  , 

I  -v  =  a(kx,k  :  z  =  0)-^—  =  jkza{kx,ky :  z  =  0)e }K\  (1.33) 

oz  oz 

Equating  (1.32)  and  (1.33)  and  evaluating  in  the  transducer  plane  at  z  =  0  yields 


jka(kx,ky :  z  =  0)  =  jcop,  J  J 


-dO0(x,y,z) 


-j(kxx+kvy) 


dxdy.  (1.34) 


Since  the  fluid  particle  velocity  vector  v(r,t)  at  any  point  in  the  background  medium  and  the  velocity 
potential  ® 0(r,t)  at  the  same  position  and  time  are  related  according  to 


v(r,0  =  -V<D(r,0, 


(1.35) 


we  can  recognize  the  factor 


d®0  (x,y,z) 
dz 


(1.36) 


appearing  in  the  integrand  of  (1.34)  as  the  component  of  the  fluid  particle  velocity  normal  to  the 
transducer  surface  in  the  plane  z  =  0.  In  other  words,  using  ( 1 .25)  and  ( 1 .36),  we  have 


vz(x,y,z  =  0)  = 


-0®o(  x,y,z) 


(1.37) 


whence  ( 1 .34)  becomes 


jk_a(kx,ky :  z  =  0)  =  jcop0  j  jv_(x,y,z  =  0)e  JiKx+kyy)dxdy.  (1.38) 
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Solving  (1.38)  for  a(kx,ky  :  z  =  0)  ,  we  have 


cc(kx , ky :  z  =  0)  =  J  J  vz  (x,  y,  z  =  0)e  Ak'x+k,y)dxdy,  (1.39) 


which  expresses  the  spatial  frequency  spectral  density  function  for  the  incident  acoustic  field  in  the 
plane  z  =  0  in  terms  of  the  component  of  the  fluid  velocity  normal  to  the  transducer  surface  in  that  same 
plane.  Insertion  of  the  spatial  frequency  spectral  density  function  of  (1.39)  into  (1.21), 


1  CO  00 

«0  (*>  y> z )  =  “  2  j  j  a(K » K  :  z  =  0 )eJ(kxX+kyy+K=)dkdky ,  ( 1 .40) 

■  '  *  —co—oo 

gives  the  acoustic  pressure  field  u0(x,y,z)  everywhere. 

For  the  transducer  surface  velocity  waveform  of  (1.25),  the  spatial  frequency  spectral  density 
function  for  the  incident  acoustic  field  in  the  plane  z  =  0  is  according  to  (1.39), 


,  .  .  .  CO  CO 

a(k  ,k  :  z  =  0)  =  — —  [  [  rect 

rect 

M 

e,il-”,’,)dxdy.  (1.41) 
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We  can  immediately  recognize  the  integral  in  (1.41)  as  the  Fourier  transfonn  of  the  two-dimensional 
rectangle  function  of  dimensions  \a,a\  in  x  and  y.  The  result  of  this  transformation  is  well  known  to 
give  a  product  of  sine  functions  thus: 


a(kx,ky  :z  =  0) 


cop0v0a 
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(1.42) 


where  sinc(x)  =  sin(x)/ x.  Inserting  (1.42)  into  (1.40)  gives  for  the  incident  acoustic  field 


cop0v0o2  }  °r  1  .  (k,a^ 

u0(x,y,z)  =  0  °2  j  J— smc 

—00—00'^'^ 


^  k,a  ^ 


smc 


V  ^  J 


V  2  J 


eKk*x+k>y+Kz)dkdkv.  (1.43) 


It  is  important  to  realize  that  the  factors  in  (1.43)  involving  k_  cannot  be  pulled  out  of  the  integrand 


since  kT  depends  on  kx  and  k  according  to  (1.19).  We  should  therefore  write  the  preceding 


expression  more  explicitly  as 
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u0(x,y,z )  = 


(0Povo 


\k20-k2x-k2y 


ka')  .  ( k  a 


arsine  ::i—  sine  — —  eJ(kxX+k,y)dkxdk  .(1.44) 


Equation  ( 1 .44)  is  the  non-unitary  inverse  Fourier  transform  of  the  bracketed  function  multiplied  by  the 
pre-factor  (op0v0 .  Thus,  we  write  ( 1 .44)  more  succinctly  as 


u0(x,y,z)  =  op0v0F 


-i  e 


ej^~k^’z  2.  (ka\.  ka  V 

=a  sine  sine  — —  >,  (1.45) 
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where  F  1  {  }  denotes  the  two-dimensional  inverse  Fourier  transfonn  with  respect  to  k  and  k  of  the 
function  in  brackets.  Invoking  the  convolution-multiplication  theorem,  we  can  write  (1.45)  equivalently 


u0(x,y,z)  =  cop0v 0  F 


^  vj  ,  2  •  (ka\.  kya\\ 

,  =  >  *  *F  <  a  sine  — 1 —  sine  — —  >  ,  (1.46) 

IK-K-K  1  2  J  l  2  JJ 


where  *  *  denotes  two-dimensional  convolution  with  respect  to  the  spatial  variables  x  and  y.  T  he 
inverse  Fourier  transform  of  the  product  of  sine  functions  is  simply  the  apodization  function  (1.26)  of 
the  transducer  velocity  waveform: 


zz-i  2  •  f  ka  \  .  f  k  a 
F  {a  sine  sine  — — 
2  2 


x  )  y 

=  rect  —  rect  — 


(1.47) 


Examining  the  first  term  in  brackets  in  (1.46),  we  see  that 


lk0  -  kx  -  k; 


j(kxx+k  y) 


kl-kl-kl 


dkdk. 
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(1.48) 


Combining  the  exponential  factors  in  ( 1 .48),  we  have 


Ik 0  -  K  -  k; 


IK-K-K 


=-  dk  dk  . 
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(1.49) 
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As  any  spherical  wave  can  be  decomposed  in  a  plane  wave  basis,  equation  ( 1 .49)  represents  the  plane 
wave  decomposition  of  a  spherical  wave  centered  at  the  origin  and  characterized  by  a  wave  vector  of 
magnitude  k0 .  In  particular,  we  have  for  z  >  0, 


1  fg7V" 
j2n  |_  r  _ 


(1.50) 


where  r  =  yx2  +  y 2  +  z2  is  the  Euclidean  distance  from  the  origin  to  the  point  (x,  y,  z).  Combing  ( 1 .46), 
(1.47),  and  (1.50),  we  have  for  the  incident  acoustic  pressure  field 


u0(x,y,z) 


~jPovof 


ejk°r 

- *  *rect 
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rect 


(1.51) 


where  we  have  used  the  relationship  co  =  Ijif  between  the  angular  frequency  of  the  transducer  vibration 
and  the  corresponding  linear  frequency.  Expression  of  u0(x,y,z)  in  tenns  of  the  linear  frequency  /is 
advantageous  since  ultrasound  transducers  are  generally  always  referred  to  in  terms  of  linear  frequency 
in  MHz.  With  the  exception  of  the  velocity  amplitude  V0 ,  all  quantities  appearing  in  (1.51)  can  easily  be 

detennined  from  knowledge  of  the  acoustic  background  medium  and  the  transducer  employed  for 
measurements. 

Equation  (1.51)  represents  the  sought-after  expression  for  the  spatial  dependence  of  the  incident 
acoustic  pressure  field  output  by  a  monochromatic  transducer  plate  of  dimensions  \a,a\  centered  on  the 
z-axis  and  lying  entirely  in  the  plane  z  =  0.  It  suggests  that  the  incident  field  is  given  quite  simply  by  the 
convolution  of  a  spherical  wave  with  a  two-dimensional  indicator  function  describing  the  transducer 
aperture.  In  reality,  (1.51)  is  simply  a  restatement  of  Huygens’  Principle,  which  considers  that  each 
infinitesimal  transducer  surface  element  acts  as  an  acoustic  monopole  emitting  a  spherical  wave.  This 
can  be  seen  by  writing  the  convolution  of  (1.51)  explicitly: 
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Since  the  rect  function  is  nonzero  only  when  the  absolute  value  of  its  argument  is  <1/2,  we  can  write 
(1.52)  as 
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Because  yj(x  —  ^)2  +(y-i])2  +  z°  is  the  Euclidean  distance  between  an  arbitrary  point  (<g, //,())  on  the 

planar  transducer  surface  and  a  field  point  of  interest  (x,  y,  z),  the  surface  integral  in  ( 1.53)  represents  the 
summation  at  the  field  point  of  an  infinite  number  of  spherical  waves,  each  of  which  is  centered  on  a 
different  infinitesimal  surface  element  of  the  transducer.  We  note  that  the  integration  limits  in  (1.53)  run 
from  \—  a  1 2,a  1 2\  in  both  the  d,  and  rj  directions  (i.e.,  the  integration  is  performed  over  the  entire 
transducer  surface).  According  to  Huygens’  Principle,  the  total  field  output  by  a  continuous  source 
distribution  is  given  by  the  superposition  of  all  spherical  waves  emitted  by  the  infinitesimal  acoustic 
monopoles  composing  the  source.  This  superposition  leads  to  regions  of  constructive  and  destructive 
interference  between  interfering  spherical  waves  and  manifests  itself  in  observable  patterns  of 
diffraction.  Ultimately,  the  cause  of  these  observed  interference  patterns  is  traceable  to  differences  in 
path  length  traveled  by  individual,  coherent  spherical  waves  in  propagating  from  their  respective  source 
centers  on  the  transducer  surface  to  a  field  point  of  interest. 

While  (1.51)  in  the  form  of  (1.53)  is  exactly  equal  to  the  well-known  Rayleigh-Sommerfeld 
diffraction  integral  for  a  square  transducer  plate  of  dimensions  [-a,«]  centered  at  z  =  0  and  excited  by 
the  velocity  waveform  (1.25),  the  efficiency  with  which  modern  computer  programs  can  perfonn 
convolution  operations  makes  the  insight  provided  by  the  angular  spectrum  derivation  well  worth  the 
effort.  Figure  2  below  shows  the  normalized  in-plane  incident  acoustic  intensity  distribution  for  a  3”  by 


3”  square  transducer  piston  operating  at  a  CW  frequency  of  3.35  MHz.  The  pressure  field  intensity  is 
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depicted  in  the  plane  parallel  to  the  transducer  surface  at  a  distance  of  2”,  and  a  background  medium  of 
water  has  been  assumed.  At  3.35  MHz,  the  attenuation  coefficient  of  water  is  only  ~  0.02  dB/cm,  so 
water  represents  a  very  good  approximation  to  a  truly  inviscid  background  medium.  We  note  that  the 
square  modulus  of  the  complex  field  given  by  (1 .5 1)  is  plotted  in  Figure  2. 

Figure  3  below  presents  the  same  information  as  Figure  2  but  shows  the  spatial  distribution  of 
acoustic  intensity  looking  at  the  plane  of  interest  head-on.  For  a  planar  detector  positioned  2”  from  the 
transducer  surface,  one  would  expect  to  observe  this  type  of  diffraction  pattern  for  an  image  acquired 
with  no  scattering  object  present. 


Y  Coordinal*  (cm) 


X  Coordinate  (cm) 
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Figure  2:  The  normalized  in-plane  acoustic  intensity  distribution  for  a  3”  by  3”  square  transducer  piston  operating  at  a  CW 
frequency  of  3.35  MHz.  The  pressure  field  intensity  is  depicted  in  the  plane  parallel  to  the  transducer  surface  at  a  distance  of 
2”,  and  a  background  medium  of  water  has  been  assumed.  We  note  that  the  square  modulus  of  the  acoustic  pressure  field 
calculated  in  (1.51)  is  plotted. 
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Figure  3:  This  figure  presents  the  same  information  as  Figure  1  but  shows  the  spatial  distribution  of  acoustic  intensity  looking 
at  the  plane  of  interest  head-on.  For  a  planar  detector  positioned  2”  from  the  transducer  surface,  one  would  expect  to  observe 
this  type  of  diffraction  pattern  for  an  image  acquired  with  no  scattering  object  present. 

We  observe  finally  that  the  complete  incident  acoustic  field,  including  both  time  and  spatial 


dependence,  is  given  by  the  product  of  (1.51)  with  the  time-harmonic  factor  e  ja>>  thus: 

u0(r,t)  =  uQ(r)e~J0,t  =  -jp0v0f 
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(1.54) 


Due  to  the  linearity  of  the  acoustic  field,  for  any  transducer  array  composed  of  multiple  square 
transducers  each  of  dimension  \a,a\  and  vibrating  in  unison  at  angular  frequency  CO,  the  total  field 
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produced  by  such  an  array  is  given  quite  simply  by  appropriately  shifting  and  adding  the  individual 
transducer  outputs  as  described  by  (1.54). 


IV.  Solution  for  the  Scattered  Field  u  (r) 


The  equation  for  the  scattered  field  component  us  (r)  from  (1.10), 


(V2  +  kl)us(r )  =  -u{r)o{r),  (1.55) 

is  an  inhomogeneous  Helmholtz  equation  with  a  source  term  on  the  right-hand  side  due  to  the  presence 
of  the  complex  object  function  o(r )  .  Thus,  as  expected,  the  scattered  field  u  .  (r)  arises  because  of  the 
perturbing  influence  of  the  scattering  object  immersed  in  the  background  medium.  Equation  (1.55) 
cannot  be  solved  directly  for  the  scattered  field,  but  a  solution  can  be  written  in  tenns  of  the  so-called 
Green’s  function  for  the  above  differential  equation.  In  particular,  the  Green’s  function,  which  we 

denote  by  g(r  |r') ,  is  a  solution  to  the  related  differential  equation 

(V2  +  K)g(r\r')  =  -8{r  -  f '),  (1.56) 

in  which  the  source  tenn  u(f)o(f )  of  (1.55)  has  been  replaced  by  the  Dirac  delta  function  8(f  —  r')  . 
Physically,  the  function  8(r  —  r')  represents  a  single  point  inhomogeneity  located  at  position  r',  and 


therefore  the  Green’s  function  g(f  r')  represents  the  scattered  wave  field  at  r  due  to  the  point 


inhomogeneity  at  r'.  In  this  sense,  the  Green’s  function  is  analogous  to  the  impulse  response  function 
of  linear  systems  theory.  F  or  a  background  medium  of  infinite  extent,  the  appropriate  three-space 
Green’s  function  satisfying  (1.56)  is  given  by 


g(r\r') 


(1.57) 
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Since  (1.57)  is  a  function  only  of  the  position  vector  difference  r—r',  we  write  the  functional 
dependence  of  g(r\r')  more  explicitly  as 

eJk olM 

-r')  =  — : - 7  *  (1.58) 

4tt|F  -  r '| 

Equation  (1.58)  corresponds  to  the  outgoing  scattered  field  generated  by  a  hannonic  point  source  at  r' 
that  radiates  a  spherical  wave;  the  wave  amplitude  falls  off  as  the  inverse  of  the  distance  between  the 
point  source  at  r '  and  the  field  point  r  of  interest. 

It  should  be  emphasized  that  the  Green’s  function  presented  in  (1.58)  corresponds  to  the 
unbounded  case  in  which  an  acoustic  monopole  radiates  a  scattered  spherical  wave  into  an  infinite 
background  medium.  In  practice,  the  background  medium  is  confined  by  rigid  barriers  that  are  capable 
of  reflecting  acoustic  waves  back  through  the  medium.  Thus,  in  the  case  of  a  bounded  medium,  the  free 
space  Green’s  function  (1.58)  may  not  provide  an  accurate  solution  for  the  scattered  wave  field,  and  a 
different  form  of  the  Green’s  function  may  be  needed  to  account  for  reflections.  Since  the  amplitude  of 
the  scattered  wave  field  in  (1.58)  varies  inversely  with  distance  from  the  scattering  center,  one  could 
expect  the  free  space  Green’s  function  to  provide  an  accurate  representation  of  the  scattered  field  when 
the  scattering  object  lies  far  from  the  boundaries  of  the  background  medium  (/.<?.,  for  weak  reflections). 

The  utility  of  the  Green’s  function  (1.58)  in  solving  (1.55)  for  an  arbitrary  object  function  o(f ) 
derives  from  the  linearity  of  the  Helmholtz  equation.  To  see  this,  we  observe  that  we  can  write  the 
source  tenn  u(r)o(r )  of  (1.55)  as  an  integral  over  delta  function  impulses,  where  each  impulse  is 
positioned  at  a  different  coordinate  r '  in  the  acoustic  medium: 

u(f)o(f )  =  J o(f')u(r')S(f  -  r')dr\  (1-59) 

Vo 

Strictly  speaking,  the  integration  in  (1.59)  is  perfonned  over  all  three-dimensional  space,  but  since  the 
object  function  is  nonzero  only  within  the  scattering  object,  (1.59)  is  also  equal  to  the  integral  over  just 
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the  physical  extent  V0  of  the  scattering  object.  Each  impulse  8(r  —  r')  in  the  preceding  expression  is 
weighted  by  the  product  of  the  total  acoustic  field  u(r ')  =  M0(r')  +  us(f ')  and  the  object  function 
o(r')  present  at  the  location  r'  of  the  impulse.  This  makes  sense  intuitively,  since  both  w(r')  and 
o(r ')  are  related  to  the  overall  scattering  power  of  the  impulse  at  r ' .  For  example,  if  o(r ')  =  0 ,  we  are 
in  the  background  medium  where  scattering  of  the  acoustic  wave  field  does  not  occur.  On  the  other 
hand,  in  regions  where  w(r')  =  0,  there  is  no  a  coustic  field  present  to  undergo  scattering,  so  the 

scattering  power  of  the  impulses  present  in  these  regions  is  zero.  We  note  that  since  o(r')  is  nonzero 

only  for  points  lying  within  the  physical  dimensions  of  the  scattering  object,  equation  (1.59)  constitutes 
a  representation  of  the  scattering  source  distribution  as  an  array  of  point  scatterers. 

The  Green’s  function  (1.58)  represents  the  solution  for  the  scattered  field  Ms(r)  due  to  an  object 
function  o(f )  consisting  of  a  single  delta  function  impulse  at  r ' .  Because  the  Helmholtz  equation  is 
linear,  we  can  solve  (1.55)  by  decomposing  the  source  tenn  u(r)o(r )  into  an  array  of  weighted 
impulses  as  in  (1.59),  solving  (1.55)  for  each  of  these  weighted  impulses,  and  adding  these  individual 
solutions  together  to  obtain  the  total  solution  to  (1.55)  for  an  arbitrary  object  function  o(f ) .  In  other 

words,  since  the  scattered  field  due  to  the  weighted  impulse  o(r')u(r')S(r  —  r')  at  r'  is  simply  the 
scaled  Green’s  function 

jk0\r-r'\ 

o(r')u(r')g(r  -  r')  =  o(f')u(f')  6  — — ,  (1.60) 

4n\ r-r'| 

we  find  that  the  total  scattered  field  at  r  due  to  the  array  of  weighted  impulses  given  by  (1.59)  is  simply 
the  superposition 

Ms(r)  =  Jo(r')M(r')g(r  -  r')dr\  (1-61) 

Vo 
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Equation  (1.61)  represents  the  solution  to  the  inhomogeneous  Helmholtz  equation  (1.55)  for  the 
scattered  field  due  to  an  arbitrary  object  function  o(f).  We  recognize  this  expression  as  a  three- 

dimensional  convolution  of  the  source  function  u(r)o(r )  with  the  Green’s  function  gir^') .  Equation 

(1.61)  physically  represents  the  total  scattered  pressure  field  at  position  r  obtained  through 
superposition  of  the  separate  scattered  fields  due  to  all  of  the  individual  point  sources  composing  the 
scattering  source  distribution.  It  should  be  remembered  that  the  acoustic  field  appearing  in  the  integrand 
of  (E61)  is  the  total  field  ll(f ')  =  M0(r')  +  Ws  if')  .  This  makes  sense  physically;  the  amplitude  of  the 

field  scattered  by  a  single  point  source  comprising  the  object  depends  on  the  total  acoustic  field  incident 
on  that  particular  point  source.  The  field  incident  on  a  given  point  source  in  the  object  consists  of  a 
contribution  from  the  incident  field  «0(r')  output  by  the  transducer  and  a  scattered  field  contribution 

us  (r ')  due  to  the  fields  scattered  by  the  other  point  sources  in  the  object. 

Equation  (1.61)  for  It  if )  is  actually  an  implicit  solution  to  the  inhomogeneous  Helmholtz 
equation,  since  the  integrand  itself  depends  on  the  scattered  field  solution  us(f): 

us(r)  =  |o(r')[M0(r')  +  «  (r')k(r  “  r')dr\  (1.62) 

Vo 

In  the  presence  of  weak  scattering,  the  real  and  imaginary  parts  of  the  incident  field  amplitude  M0(r') 
are  much  larger  than  the  corresponding  parts  of  the  scattered  field  amplitude  us(r')  for  all  r' .  This 
allows  us  to  neglect  the  scattered  field  component  of  the  total  field  in  (1.62)  to  obtain  an  explicit 
solution  for  Usif)  in  tenns  of  the  complex  object  function,  the  incident  acoustic  field  (/.<?.,  the  field 
calculated  in  (1.51)),  and  the  Green’s  function  (1.58): 

us(r)  =  uB(r )  =  jo(r')w0(r')g(r  -  r')dr'.  (1.63) 

Vo 


39 


Equation  (1.63)  expresses  the  so-called  first-order  Bom  approximation  for  the  scattered  field  component 
of  the  total  acoustic  field.  The  Bom  approximation  is  iterative  in  nature,  as  substitution  of  (1.63)  for 
Ms(r)  back  into  (1.62)  can  improve  the  accuracy  of  the  initial  estimate  for  u  (r) .  This  procedure 
would  yield  the  second-order  Bom  approximation  for  the  scattered  field,  which  itself  could  be 
substituted  back  into  (1.62)  for  Us  (r)  to  obtain  yet  another  improvement  in  the  estimate  of  the  scattered 
field.  In  principle,  this  process  can  be  repeated  an  arbitrary  number  of  times  to  obtain  the  z^-order  Bom 
estimate  for  the  scattered  field  amplitude  Us(r ) .  The  ith- order  Bom  approximation  for  the  scattered 
field  is  given  by 

uj'\r)  =  \°(F){u0(F)  +  uBu  l) (r'))g(r  -r')dr' ,  (1.64) 

where  UB(°\r')  =  0.  Equation  (1.64)  is  identical  to  (1.62),  with  the  exception  that  Us(f ')  has  been 

replaced  by  the  estimate  UB  ’ (f7' ) .  For  our  purposes,  particularly  those  related  to  tomographic 

reconstmction,  we  will  find  the  first-order  Born  field  as  given  by  (1.63)  to  be  the  most  useful  of  the 
iterative  Born  approximations. 

Before  considering  how  the  first-order  Born  approximation  can  be  used  for  tomographic 
reconstmction,  it  is  instmctive  to  examine  (1.63)  more  closely  to  understand  the  fundamental  physical 
assumptions  underlying  the  Bom  theory.  Based  on  our  previous  discussion  concerning  the 
decomposition  of  a  continuous  scattering  distribution  into  an  array  of  point  scattering  centers,  we  see 
that  (1.63)  suggests  that  the  field  incident  on  each  point  scatterer  in  the  object  is  simply  w0(r'),  the 

incident  acoustic  field  output  by  the  transducer.  Thus,  (1.63)  implies  that  the  incident  field  output  by  the 
transducer  remains  essentially  unchanged  as  it  passes  through  the  scattering  object.  For  a  homogeneous 
cylinder  of  radius  a  having  an  acoustic  refractive  index  deviation  ns(r )  =  ns  relative  to  the  refractive 


40 


index  of  the  background  medium,  it  is  possible  to  show  that  the  Born  assumption  of  a  negligible 
perturbation  of  the  incident  transducer  field  is  satisfied  under  the  condition 


/l 

anx  <  — , 

*  4 


(1.65) 


where  A,  denotes  the  wavelength  of  the  acoustic  wave  field  in  the  background  medium.  Therefore,  a 
distinctive  feature  of  the  first-order  Born  approximation  is  that  the  conditions  for  its  validity  depend  not 
only  on  the  refractive  index  of  the  imaged  object  but  also  on  its  size. 

Because  of  its  conduciveness  to  tomographic  reconstruction,  from  this  point  forward,  we  will 
consider  only  the  first-order  Bom  approximation  for  the  scattered  field, 

us(r)  =  uB(r)=  \o(r')u0(r')g(r  -r')dr\  (1.66) 


where  the  incident  field  U0(f ')  is  the  output  of  the  acoustic  transducer  that  we  calculated  in  (1.51), 
namely 


u0(x',y',z')  =  -jp0vj 


?jKr' 


*  *rect 


rect 
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(1.67) 


In  a  single  transmission  image  of  an  object,  ultimately  what  is  measured  in  the  detector  plane  is  the 
intensity  of  the  total  acoustic  field  u(f)  ,  which  includes  components  due  to  both  the  transducer  output 
as  well  as  scattering  from  the  object  being  imaged: 

u(r)  =  u0(r)  +  us(r).  (1.68) 

By  adding  (1.66)  and  (1.67),  we  have  a  predictive  model  for  the  spatial  variation  in  acoustic  intensity 
1(f)  that  would  be  measured  in  the  detector  plane  for  an  arbitrary  object  being  imaged: 
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(1.69) 
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V.  Algorithms  for  Ultrasound  Tomographic  Reconstruction 

The  goal  of  any  tomography  experiment  involving  the  imaging  of  a  three-dimensional  object  is 
accurate  volumetric  reconstruction  of  the  object  from  a  series  of  two-dimensional  planar  projections.  In 
ultrasound  diffraction  tomography,  since  spatial  variations  in  the  acoustic  refractive  index  are 
responsible  for  wave  field  scattering,  what  we  are  actually  reconstructing  is  a  three-dimensional 
mapping  of  the  real  and  imaginary  parts  of  the  acoustic  refractive  index  of  the  object  immersed  in  the 
background  medium. 

Fundamental  to  diffraction  tomography  is  the  Fourier  Diffraction  Theorem  (FDT),  which  is 
analogous  to  the  Central  Slice  Theorem  of  x-ray  computed  tomography.  In  fact,  by  use  of  fractional 
Fourier  transforms,  it  can  be  shown  that  both  theorems  are  equivalent  in  the  short-wavelength  limit.  The 
FDT  relates  the  two-dimensional  Fourier  transfonn  of  the  scattered  component  Ms(r)  of  the  total 
acoustic  field  measured  in  the  detector  plane  to  the  three-dimensional  Fourier  transform  of  the  object 
function  itself.  The  FDT  is  extremely  powerful  and  enables  the  Fourier  space  of  the  object  function  to  be 
built  up  from  a  series  of  two-dimensional  planar  projections.  Fourier  inversion  of  this  result  provides  an 
estimate  of  the  object  function  directly.  It  should  be  realized  that  the  FDT  is  valid  only  in  the  presence  of 
weak  scattering  by  inhomogeneities  in  the  imaged  object. 

The  measurement  geometry  of  ultrasound  diffraction  tomography  consists  of  an  acoustic  source 

and  a  planar  detector  that  are  diametrically  opposed  and  positioned  on  oppos  ite  sides  of  the 

inhomogeneous  object  being  imaged.  For  ultrasound  illumination  of  the  object  o(r )  using  a  single 

plane  wave,  the  FDT  states  that  the  two-dimensional  Fourier  transform  of  the  forward  scattered  field 

Ms(r)  measured  over  the  detector  plane  is  equal  to  a  hemispherical  surface  in  the  three-dimensional 

Fourier  space  of  the  object  function  o(f )  .  This  surface  passes  through  the  origin  of  the  object  Fourier 

space  and  is  oriented  opposite  the  propagation  direction  of  the  probing  plane  wave  field.  Moreover,  the 
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radius  of  the  hemispherical  surface  is  determined  by  the  wave  vector  magnitude  (i.e.,  the  wavelength)  in 
the  background  medium  of  the  ultrasound  energy  used  to  illuminate  the  object.  By  rotating  the  detector 
and  source  assembly  around  the  object  and  imaging  at  different  angular  positions,  the  Fourier  space  of 
the  object  can  be  built  up  using  the  hemispherical  surfaces  provided  by  the  distinct  source-detector 
positions.  In  practice,  a  number  of  discrete  angular  positions  of  the  source  and  detector  are  employed, 
and  appropriate  interpolation  techniques  are  then  used  to  estimate  the  values  of  the  Fourier  transform  of 
the  object  function  where  measurement  data  is  not  available.  The  interpolated  Fourier  space 
representation  of  the  object  is  then  inverted  to  obtain  an  estimate  of  the  object  function  o(f )  . 

For  an  arbitrary  source  of  acoustic  waves,  the  FDT  as  outlined  in  the  preceding  paragraph  is  not 
directly  applicable,  since  the  incident  pressure  field  may  not  be  a  simple  plane  wave.  In  particular,  the 
incident  acoustic  wave  field  relevant  to  our  ultrasound  system  is  given  by  the  convolution  of  (1.51).  We 
recall,  however,  that  the  expression  we  derived  in  (1.51)  for  the  incident  field  was  based  on  an  angular 
spectrum  decomposition  method  involving  a  superposition  of  plane  wave  solutions  to  the  homogeneous 
Helmholtz  equation.  Thus,  it  is  reasonable  to  expect  due  to  the  linearity  of  the  Helmholtz  equation  that 
alternative  methods  based  on  the  FDT  exist  for  filling  the  Fourier  space  of  an  object  using  arbitrary 
wave  fields.  Such  algorithms  have  previously  been  developed,  and  we  present  one  of  these  methods  in 
this  section.  The  advantage  of  this  procedure  is  that  it  requires  only  two  angular  orientations  of  the 
object  separated  by  90°.  It  is  based  on  acquiring  measurements  of  the  scattered  field  throughout  the 
detector  plane  as  the  position  of  the  acoustic  source  is  varied  in  the  transducer  plane.  Obtaining 
diffracted  projections  of  the  object  in  this  manner  has  the  added  benefit  of  more  closely  resembling  the 
B-mode  scanning  procedure  of  conventional  pulse-echo  ultrasound  imaging,  a  procedure  with  which 
most  everyone  is  familiar. 

The  development  of  the  aforementioned  reconstruction  method  involves  writing  the  volume 
integral  (1.63)  representing  the  first-order  Born  approximation  in  a  way  that  relates  the  Fourier 
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transform  of  the  object  function  to  the  Fourier  transform  of  the  scattered  field.  This  procedure  can  be 
perfonned  by  making  use  of  the  angular  spectrum  method  outlined  in  section  II.  In  particular,  we 
showed  previously  that  for  an  ultrasound  transducer  centered  at  the  origin  of  coordinates  in  the 
transducer  plane,  the  angular  spectrum  decomposition  of  the  incident  acoustic  field  U0(r )  in  any  plane 
of  constant  z  is  given  by 

1  OO  CO 

u0(r)  =  —  J  J a(kx,ky :  z  =  0)eJ>k'x+kyy+Kz,dkxdky,  (1.70) 

■  '  ' *  —  CO—  OO 

where  all  variables  are  defined  as  they  were  in  section  II  and  we  recall  that  k_  —  ^ jk 02  —  kx  —  k  ]  .  If  we 

consider  a  translation  of  the  transducer  surface  from  the  coordinate  origin  of  the  transducer  plane  to  the 
point  (xt,yt)  in  that  plane,  then  it  follows  from  (1.70)  that  the  angular  spectrum  decomposition  of  the 

transducer  field  u0(r  \xt,yt)  is  given  by 


«0 (F :  xi’yt) =  y  j  j , ky :  z  =  0)e 

■  *' *  —CO— CO 


j(kx  (x-x,  )+k  (y-y,  )+k.z) 


dkxdkv .  (1.71) 


Equation  (1.71)  represents  a  simple  shift  of  the  entire  incident  field  expressed  in  (1.70)  by  Xt  in  the  x- 


direction  and  yt  in  the  v-direction.  This  is  just  what  we  expect,  since  the  only  physical  difference 
between  the  situations  described  by  (1.70)  and  (1.71)  is  a  relative  shift  (xt,yt)  of  the  transducer 
position  in  the  source  plane. 

Considering  now  the  free  space  Green’s  function  given  by  (1.58),  we  observed  earlier  that  a 
spherical  wave  centered  at  the  origin  and  characterized  by  a  w  ave  vector  of  magnitude  k0  has  the 
angular  spectrum  decomposition  shown  in  (1.50),  namely 
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where  z  >  0.  From  (1.72),  it  follows  that  the  angular  spectrum  decomposition  of  the  Green’s  function 
(1.58)  is  given  by 
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(1.73) 


where  k0  =  (a,j3,^jk0  —  a  —  (3  ) . 

If  we  insert  the  decompositions  (1.71)  and  (1.73)  into  the  first-order  Born  approximation  (1.63), 
we  obtain  an  approximate  expression  for  the  scattered  field  when  the  transducer  is  centered  at  (x, ,  yt ) 
in  the  transducer  plane: 


us(r:xt,yt)  = 
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x  (1.74) 


j(kx (x'-x, )+k  (y'-y, )+kzz') 


dkxdky , 


dr'. 


For  a  p  lanar  detector  lying  in  the  plane  z  =  D,  the  scattered  field  amplitude  at  the  position 
r  =  {xr,yr,D)  in  the  detector  plane  is  therefore  given  by 


us(r:x„yt)  = 
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x  (1.75) 


\ \a(kx,ky:z'=0)e 


j(kx (x'-x, )+k  (y'-y,  )+kxz') 


dkxdky , 


dr'. 


Letting  xt  — >  —  xr  and  yt  — >  —yi  in  (1.75)  gives 


K(f:~xt-yt)  = 
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Taking  the  two-dimensional  Fourier  transfonn  of  both  sides  of  (1.76)  with  respect  to  the  spatial 
coordinates  (xt,yt)  of  the  transducer  position,  and  then  taking  the  two-dimensional  Fourier  transform 
of  both  sides  of  the  resulting  equation  with  respect  to  the  spatial  coordinates  ( Xr,y  r )  of  the  field  point 
in  the  detector  plane  gives 


U3(a,P\kx,ky)  = 
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where  Us(a,  x,ky)  of  (1.76)  is  related  to  the  scattered  field  amplitude  Us(f  :—xt,—yt)  in  the 
plane  z  =  D  according  to 


u, {ct,fd ;kx, ky )  =  J  J  J  J us (xr ,yr,D:-x/ ,-yt )e 
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x  (1.77) 
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Recognizing  the  integral  in  (1.76)  as  the  three-dimensional  Fourier  transform  of  the  object  function 
o(r')  evaluated  at  (ct  —  kx,j3  —  kv,y  —  k_ )  enables  us  to  write  (1.76)  as 


Us(a,frkx,ky) 


— a(kx ,  k  :  z'  =  0 )0(a  -  kx ,  p  -  k  ,  y  -  kz ), 
2  y 


(1.78) 


where  y  =  J k(j  —  a  —  (3  and  k,  —  Jk0  —  k~  —  k~ .  It  can  be  shown  using  (1.78)  that  the  scattered 


field  data  obtained  for  a  single  projection  of  the  object  yields  coverage  in  the  frequency  domain  of  the 
object  function  defined  by  two  spheres.  A  single  projection  view  in  this  imaging  procedure  corresponds 
to  a  single  orientation  of  the  object  relative  to  the  source  plane  but  multiple  transducer  positions  in  the 
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source  plane.  For  a  single  projection  view,  the  scattered  field  is  collected  in  the  detector  plane  as  a 
function  of  the  transducer  position  in  the  source  plane.  By  rotating  the  object  90°  and  performing  the 
same  scanning  procedure  of  the  transducer  in  the  source  plane,  it  is  possible  to  generate  the 
complementary  spheres  in  the  frequency  domain  of  the  object  function  and  thus  fill  in  the  Fourier  space 
of  the  object. 


47 


